Instabilities of interacting electrons on the honeycomb bilayer 
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We investigate the instabilities of interacting electrons on the honeycomb bilayer by means of 
the functional renormalization group for a range of interactions up to the third-nearest neighbor. 
Besides a novel instability toward a gapless charge-density wave we find that using interaction 
parameters as determined by ab-initio calculations for graphene and graphite puts the system close 
to the boundary between antiferromagnetic and quantum spin Hall instabilities. Importantly, the 
energy scales for these instabilities are large such that imperfections and deviations from the basic 
model are expected to play a major role in real bilayer graphene, where interaction effects seem to 
be seen only at smaller scales. We therefore analyze how reducing the critical scale and small doping 
of the layers affect the instabilities. 

PACS numbers: 71.10.-w, 73.22.Pr 
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Recently, a number of experiments have addressed the 
ground state properties of bilayer graphenefll-Q, trying 
to clarify possible signatures for the interaction between 
the electrons. So far, however, the precise nature of the 
ground state and whether the spectrum is gapped or gap- 
less remain under controversial discussions. A number of 
competing instabilities of the semi-metal, e.g. toward a 
gapped antiferromagnetic (AF) state, a gapless nematic 
phase and even toward gapped topological phases as the 
quantum anomalous and quantum spin hall states are de- 
bated and compared to the experimental findings. These 
proposals are subject of numerous theoretical works 
14l |. Comprehensive renormalization group (RG) investi- 
gations of the different interaction processes giving rise 



part of the Hamiltonian within layer I has the form[lf 



13| 



to the various instabilities have been performed [12| 
within a continuum model for the vicinity of the band 
crossing points (BCPs), also addressing the role of the 
range of the interactions. These studies indicate that the 
ground state properties may depend decisively on the 
profile of the interaction. In this context it is impor- 
tant to notice that the effective interaction parameters 
and their spatial dependence for the usual low-energy 
models of graphene and graphite have been calculated 
by ab-initio techniques [l5|. using the constrained random 
phase approximation (cRPA) that takes into account the 
screening due to bands further away from the Fermi level. 
In this work we use a functional RG (fRG) scheme for an 
unbiased investigation of the instabilities of the A-B or 
(Bernal-)stacked bilayer honeycomb lattice. Besides ex- 
ploring the general picture, we use the interacti on p aram- 
etcrs as specified by the ab-initio calculations [15[. This 
allows us to get closer toward a realistic picture of the 
dominant ordering tendencies in bilayer graphene. 

We start with the model for the graphene bilayer at the 
charge neutrality point. With respect to a site in sublat- 
tice A in honeycomb layer I, the three nearest neighbor 
(n.n.) sites on sublattice B in the same layer are given 
by the shifts 5i with t G {1,2,3}. The n.n. hopping 



^ 5i)ai^s{R)+\\.c^ , where 
,{R) the annihilation op- 



s —^,\. is the electron spin, a;. 

erator of an electron at site R on sublattice A, b\ ^ the 
creation operator on sublattice B, and I the layer index. 
The relative sign for the neighbors 5i in the two layers is 
due to the A-B-stacking. The interlayer hopping t±^ be- 
tween sites that lie on top of each other (Ai and A2-sites) 
reads0 = i-L i? {a\,s{R)a2,s{R) + h.c). Diago- 



nalizing iJfrcc = E; 



=1,2^; 



-\-H^ results in 4 bands. Two 



of them have two incquivalcnt BCPs K and K' at the 
Brillouin zone (BZ) corners at the Fermi level. While 
i_L = yields two copies of the single layer with linear 
dispersion (Dirac cones) at K and K' , tj_ 7^ results in a 
quadratic dispersion near those Fermi points. This allows 
interaction-driven instabilities for arbitrarily small cou- 
plings. The instabilities of a sin gle quadratic band cross- 
ing point were analyzed in Refs. IT], ll8[ . We ignore trigo- 
nal warping terms for the time being, but comment later 
on the effect of perturbations on the quadratic BCPs. 

As interactions, we account for an onsite repulsion U , 
a n.n. density-density interaction Vi and a second-n.n. 
density-density interaction V2. For these terms, cRPA 
values are listed in Ref. [l3| • In addition, for checking the 
robustness of our results, we consider a third-n.n. repul- 
sion V3. The interaction Hamiltonian reads 

H\nt = U'^n^^-^rii^i + Vi ^ nj,snj,s' (1) 

{{i-J))-s,s' {{{i,j))),s,s' 

Here i,j run over the lattice sites, but pairs are included 
only once. The unitary transformation from the orbital 
fields to the bands diagonalizing iffroc has to be carried 
out on iJint as well. This adds orbital makeup to the 
interaction terms in band representation, leading to an 
additional angular dependence of the interactions near 
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FIG. 1: Left panel: Sketch of Bernal-stacked bilayer graphene 
with the sublattices A;,B; and layer index I £ {1,2}. Right 
panel: Patching scheme of the Brillouin zone in the fRG. The 
black dots denote the momentum vectors at which the cou- 
pling function is evaluated. 



the K, if'-points. In order to resolve this dependence 
we use an angular patching of the interaction terms as 
explained below. This wa y o ur study goes beyond the 
'g-ology' approach in Ref. [l3| . 

We employ a functional renormalization group (fRG) 
approach for the one-particle-irreducible vertices of a 
fermionic many-body system (for a recent review, see 
[igj). In this scheme, an infrared regulator with energy 
scale A is introduced into the bare propagator. The RG 
flow is generated upon variation of A. By integrating the 
flow down from an initial scale Aq of the order of the 
bandwidth to the infrared A — ^ 0, one smoothly interpo- 
lates between the bare action of the system and the effec- 
tive action at low energy. Here, the hierarchy of flowing 
vertex functions is truncated after the four-point (two- 
particle interaction) vertex. This vertex is described by 
a coupling function Vj\(fci, ^2; fcs, fc4) with a discretized 
momentum dependence for incoming quantum numbers 
ki , ^2 and outgoing k^jk^. Here ki includes a wavevcctor 
ki, a Matsubara frequency w^, a spin projection Si, and 
for our case also band-, or orbital indices. For simplicity, 
and as one is interested in groundstate properties, exter- 
nal frequencies are set to zero. Furthermore, in order to 
keep the calculations doable, we neglect self-energy cor- 
rections. This approximate fRG scheme then amounts to 
an inflnite-order summation of one-loop particle-particle 
and particle-hole terms of second order in the effective in- 
teractions. It allows for an unbiased investigation of the 
competition between various correlations, by analyzing 
the components of Va(A;i, A:2; fca, A;4) that create instabil- 
ities by growing large at a critical scale Ac[l3|- With 
the approximations mentioned above, this procedure is 
wcU-controUcd for small interactions. At intermediate 
interaction strengths we still expect to obtain reliable re- 
sults. In any case, the fRG takes into account effects 
beyond mean-field and RPA. The scale Ac can be inter- 
preted as estimate for ordering temperatures, if ordering 
is allowed, or at least as temperature below which the 
dominant correlations should be clearly observable. Fur- 
thermore, one can understand Ac as energy scale for the 



modification of the spectrum, typically by a gap. 
The discretization of the interaction vertex Va is imple- 
mented by dividing the BZ into A'' patches with constant 
wavevector-dependence within one patch. The represen- 
tative momenta for the patches are chosen to lie at or 
close to the Fermi level. The patching scheme is shown 
in Fig. [1] with TV = 24. Each of the four momenta in 
Va(A:i, ^2; fca, ^4) is additionally equipped with a sublat- 
tice index and a layer index. Momentum conservation 
fixes one of the four wavevectors. Altogether this results 
in a 4^ X component coupling function Va- 
In this work we study the flow at temperature T = 0. We 
find flows to strong coupling with non-zero critical scales 
Ac for all choices of non- vanishing interaction terms pro- 
vided t± ^ 0, due to the non- vanishing density of states 
at the Fermi level of the coupled layers In practice, 
the flow is stopped at some flnite value for the largest in- 
teraction component of the order of twice the bandwidth. 
Then Ac is determined by extrapolation. 
Here we present fRG results at zero temperature for inter- 
layer hopping t± = O.li as estimate for values in the liter- 



ature, e.g. for graphite[21[ or for few-layer graphene[22[. 
We also analyzed larger t± < t, without major qualita- 
tive differences. The monolayer-case t \ = was studied 
by fRG in Rcfs. [H, and in Rcfs. [H ^ for larger 
The next-neighbor hopping t « 2.8eV sets the 



2^ and in Rcfs 

doping 

energy unit. We then study the parameter space spanned 
by U, V\ and V2 up to the cRPA parameters found in[l5|. 
By identifying the leading tendencies, i.e. the strongest 
class of divergent couplings, we encounter rich tentative 
phase diagrams shown in Fig. [2j The drawn boundaries 
arc guides to the eye. Typically, the flows change contin- 
uously from one regime to the other without drastic fea- 
tures in Ac. Hence, while without including self-energy 
effects a possible suppression of the critical scales due to 
quasiparticle degradation is not captured, we expect that 
these transitions are of first order. We now discuss the 
various ordering tendencies found for given parameters 
and how they are revealed in the fRG flow. 
Antiferromagnetic (AF) spin density wave (SOW) insta- 
bility. In the fRG data, the flow towards the AF-SDW 
is seen as a leading divergence of interaction components 
with zero momentum transfer in the spin channel. It 
features an attractive sign for intra-sublattice scatter- 
ing and a repulsive sign for inter-sublattice processes, 
in complete correspondence to the single layer [2 3|. The 
interlayer sign structure can be read off from the RG 
data as well. In detail, the leading part of effective 
interaction in this case (and ignoring frequency depen- 
dences) reads Haf — ~Tf^oo'^oo'^o£o'S° ■ S° with 

Voo' > and 5^° = ^ Es.^.s' '^-^'ct^.o^.' - ^"^ 



2 ^k,s,s' 

depend on the orbital, = +1 for o G {01,62} and 
60 = —1 for o S {02, &i}. The effective interaction has 
become infinitely-ranged due to the sharpness in momen- 
tum space. A mean-field decoupling of Haf results in an 
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AF spin alignment in each layer where a net spin (e.g. 
'up') moment is located on the Ai- and B2-sublattices, 
and an opposite net spin ('down') moment on the Bl- 
and A2-sublattices. In absence of spin-orbit interactions 
the spin quantization axis is not fixed. This phase is 
sometimes also referred to layer antiferromagnet 

(LAF) 3. A closer look at the fRG data shows that 
the intralayer components Voo' on the B sublattices grow 
faster toward the instability than the couplings on the 
A-sublattices, pointing to a larger spin moment on the 
B sublattice, in agreement with Quantum Monte Carlo 
(QMC) simulations pol for the same system with pure 
onsite interactions. 

Charge density wave (CDW) instability. Here we en- 
counter diverging interactions in the density channel, 
again with zero momentum transfer, with opposite 
signs for the intra- and interorbital interactions. In 
detail, we observe the effective interaction i?cDW = 
-wEo.o'^oo'eoeo'iV°Af°' with Voo' > and N" = 

y^r ct cr . Within a layer, this results in an 
infinitely-ranged attraction for sites on the same sublat- 
tice and repulsion for sites on different sublattices. The 
sign-structure between the layers favors an enhanced oc- 
cupancy of the Ai and B2 sublattices and a reduced oc- 
cupancy on the Bi and A2 sublattices or vice versa. The 

electronic spectrum becomes gapped by this ordering^ 

Quantum Spin Hall ( QSH) instability. The QSHj27l |28|| 
phase breaks spin-rotational symmetry, whereas time re- 
versal symmetry remains conserved. In the fRG flow, 
spin interactions with zero wavevector transfer diverge, 
with an additional sign structure that alternates between 
K and K' points, and between the sublattices. The effec- 
tive interaction reads i/qsH = - jj Eo.o' ^oo'Eoe, 
with Voo' > and S] = 5 Efe,.,,' 4'?. 
including a /-wave form factor — sm(kx) — 



QO CO 

'^f^f 

k,s,o k,s',o 



k \/3 

2sin(-^)cos( 2 )• III 9- mean-field treatment of this 
effective interaction, a purely imaginary Kane-Melepvj 
order parameter Aij = J2s.s' '^ss' {al ,,aj^s') = A* j with 
next-nearest neighbors i,j is induced. In the fRG, the 
chirality of the state comes out the same in the two lay- 
ers for the same spin, i.e. there should be two edge modes 
with the same propagation direction per spin. Hence, the 
edge state would not be topologically protected [isj. 
Three- sublattice CDW instability (CDW3). Interestingly, 
we found another instability for smaller U and V2/t < 
1.0, which so far seems not to be mentioned in the liter- 
ature. Here a site-centered CDW tendency with a finite 
momentum transfer Q ^ K — K' = K' grows during the 
fRG flow. The wavevector dependence of Va near Ac for 
this CDWa-instability is shown in Fig. [31 The sharp fea- 
tures belong to wavevector transfer ±Q. Only processes 
with initial ki and k2 near different BCPs grow strongly, 
because only then the final states after scattering by ±Q 
lie near the BCPs as well. This causes the interruption of 
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FIG. 2: Tentative £RG phase diagrams for U/t = 0,1,2,3.5 
and t± — O.lt. The area encircled by the dashed line shows 
the region of the cRPA parameters of Ref. p^ . 



the horizontal features in Fig. [3l The corresponding ef- 
fective interaction given by these leading terms near the 
divergence becomes 



HcDWs 



with 

Q 



N ^ 



eoeo'{NZN". + K.Nt)i2) 



^k. 



•■ k+Q,s,o k,s,o' 



In this equation it is 



understood that the wavevectors k range in the vicinity 
of the K and K'- points, as only there the interactions 
grow large. In a variational treatment, ^ is minimized 
by complex expectation values (N^) = eo|Ao|e*". These 
break the translational symmetry by density modulations 
cx cos{Q ■ R+a). Based on the fRG data for Voo', the |Ao| 
should be of comparable magnitude. Each original sub- 
lattice o is broken up into three sublattices (see Fig. [3|). 
Changing a reorganizes the charges within the three new 
sublattices while keeping the average constant. The dis- 
crete rotational symmetry is broken completely for gen- 
eral a. Hence this state should be observable directly in 
scanning tunneling experiments. The fermionic spectrum 
is gapless and has 4 Dirac cones near the center of the 
new reduced BZ. 

Let us now discuss the relation to bilayer-graphene. 
In Ref. [Til the interaction parameters for single-layer 
graphene and graphite were estimated by ab-initio meth- 
ods. We expect bilayer graphene to interpolate between 
these cases. The area of these ab-initio interaction pa- 
rameters is shown as a dashed line in Fig. Eljd). In the 
fRG for this parameter range SDW and QSH instabili- 
ties compete. Which one is leading depends however on 
the precise values of the interaction and also the hopping 
parameters, so that we refrain from giving a definite pre- 
diction. We have also checked that a third-n.n. repulsion 
V3 < t does not change the nature of the ground state for 
these parameters. As both states have a non-zero single- 
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FIG. 3: CDW3 instability. Left & middle plot: EfTective 
coupling in units of t versus incoming patch indices ki and 
k2 for U = 0,Vi = 0.5t,V2 = 0.5t and t± = O.lt, with first 
outgoing wavevector fca on patch 1. Left plot: 01=02 = 
03 = 04. Middle plot: oi = 03 and 02 — 04 ^ o\. Right: 
Qualitative charge distribution in the CDW3 mean-field state 
for ct — 1/6 and equal |Ao|s in one layer. The size of the 
symbols indicates the charge density, the three sizes in one 
(original) sublattice average to 1. 




n/t a 



FIG. 4: Left panel: fRG critical scale Ac versus chemical po- 
tential n for graphene (squares) and graphite (triangles) cRPA 
interaction parameters [1^. Right panel: Critical scale Ac as 
function of a rescaling parameter a for graphene (squares) 
and graphite (triangles) cRPA interaction parameters. 



particle-gap, they are compatible with the experimental 
spectrum of Ref . Q ■ A distinction might be made from 
testing for time-reversal symmetry breaking, which only 
occurs for the SDW state. As the QSH state corresponds 
to two copies of the Kane-Mele single-layer state, gapless 
edge state transport might be spoiled by impurities psf. 
Further, for the /-wave order parameter of the QSH, im- 
purities will be detrimental, and the experimental gap is 
only seen in ultraclean samples Q. 

We also investigated the impact of non-vanishing 
chemical potentials (/i 7^ 0) to mimic the effect of im- 
purities or small dopings on the groundstate, with focus 
on the cRPA interaction parameters and t± = O.lt. In 
the range between fj, € (0,0. 5t] the critical scale Ac only 
changes mildly and the groundstate remains unchanged. 
Note that the instability scales deduced for these realistic 
parameters are huge (up to t) , due to the perfect nest- 
ing between the two bands forming the quadratic BCPs. 
We cannot guarantee that our method still works quan- 
titatively in this regime. Most likely the inclusion of self- 
energy corrections and frequency dependence of the in- 
teractions will reduce the fRG scales. However, we show 
this data on purpose in order to expose clearly that the 
experimental energy gapsQ 2meV ~ 10~H are sev- 
eral orders of magnitude smaller than the gap scales one 



gets theoretically in this simple modeling, using a method 
that should be expected to be more realistic than mean- 
field theory or simpler perturbative arguments. The 
single-layer system for pure onsite interactions offers a 
possibility to compare the energy scales in the fRG with 
published non-perturbative QMG data. At U ~ 4.3t 
where the AF order sets in in QMG [2^ above a narrow 
spin- liquid regime, the single- particle gap is ~ O.lSt. The 
temperature- flow fRGfisj with upper bandwidth ±t 
gave Ac ~ 0.16t, with a critical Uc ^ 3.8t for AF order. 
Our present code with energy-cutoff and integration over 
the full bandwidth gives a higher Ac = 0.85t a.t U ~ 4.3t, 
which is already far above the corresponding Uc = 2.8t 
within this scheme. Hence the overestimate may reach 
factors ~ 5 — 6, but in a limited range of parameters near 
the opening of the gap at Uc in QMG. In the bilayer sys- 
tem, Uc is expected to be zero and high scales occur for a 
wide range of parameters. We expect that other methods 
will confirm this fRG result. Of course, in the experimen- 
tal system, unintentional doping and potential variations 
will lead to a reduction of the energy scales. But as one 
can see from the doping dependence just described, a sin- 
gle perturbation has to be rather strong to be effective, 
and one may actually need a combination of factors like 
doping, disorder, trigonal warping, additional reductions 
of the interaction parameters etc., to reduce the scales 
down to values compatible with the experiment. In or- 
der to estimate the dependence of the ground state on the 
critical scale, we introduced a global rescaling parameter 
a for the interaction terms, i.e. U aU, V\ -> aVi, 
V2 — !■ aV^- We found that a does not change the nature 
of the ground state, i.e. for the graphene parameters we 
observed a QSH instability for all a and for the graphite 
parameters the SDW instability. Ac decreases by two 
orders of magnitude as a < 0.5, as shown in the right 
plot of Fig. m The marked reduction is caused by the 
strongly energy-dependent density of states. 

In summary, we have presented a fRG study of interac- 
tion driven instabilities in the honeycomb bilayer model. 
Besides a novel gapless GDW state, we found that us- 
ing ab-initio estimates for the band-structure and non- 
local interaction parameters for bilayer graphene leads 
to a narrow competition of quantum-spin-Hall and AF- 
SDW instabilites, where details might decide what the 
actual groundstate is. Another important information is 
the typically high energy scale for the breakdown of the 
gapless state in the model used by us and many other 
theoretical studies. At present, these high scales do not 
seem to be reflected in the experiments, and more re- 
search is needed to understand this discrepancy. 
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